clc;clear
%%% G2010
load h_rms;H=z;H(H==-999)=nan;H2d=H;H=H(~isnan(H2d));
load ks;ks=z/1e3;ks=ks(~isnan(H2d));
load kn;kn=z/1e3;kn=kn(~isnan(H2d));
load nu;nu=z;mu=2*nu+2;mu=mu(~isnan(H2d));
load theta;theta0=z;theta0=theta0(~isnan(H2d));
fcori=zeros(5400,2250); % size G2010
fcori1=abs(2*2*pi/(24*3600)*sind(y));
for i=1:5400
    fcori(i,:)=fcori1;
end
fcori=fcori(~isnan(H2d));clear fcori1
P0=2*pi.*H.*H.*(mu-2)./kn./ks;
%%% wavenumber
kf=-10.^(-6:.1:-1);kf=kf(end:-1:1);kf(end+1)=0;kf(end+1:end+length(kf)-1)=10.^(-6:.1:-1);
k=0.5*(kf(2:end)+kf(1:end-1));dk=kf(2:end)-kf(1:end-1);clear kf
l=k(52:end);dl=dk(52:end);
for i=month
    N=rdda(['N_goff',num2str(year),'.bin'],[5400,2250],i,'real*4','b');
    G2d=zeros(5400,2250,30);Fr2d=G2d;a=1;
    for j=1+num(i):num(i+1)
        u=rdda(['u_goff',num2str(year),'.bin'],[5400,2250],j,'real*4','b');u=u(~isnan(H2d));
        v=rdda(['v_goff',num2str(year),'.bin'],[5400,2250],j,'real*4','b');v=v(~isnan(H2d));
        ubar=(u.^2+v.^2).^.5;
        E0=1027.5/4/pi/pi;E=zeros(length(H),1);P_prop=E;
        parfor m=1:length(k)
            for n=1:length(l)
                theta=atand(k(m)/l(n));
                klsqr=k(m).^2+l(n).^2;
                Udotk=u.*k(m)+v.*l(n);
                P=P0.*(klsqr.*((cosd(theta-theta0)).^2)./kn.^2+...
                    klsqr.*((sind(theta-theta0)).^2)./ks.^2+1).^-(mu/2);
                P(Udotk.^2>N.^2 | Udotk.^2<fcori.^2)=0;
                E=E+E0.*abs(Udotk)./klsqr.^.5.*P.*(Udotk.^2-fcori.^2).^.5.*(N.^2-Udotk.^2).^.5.*dk(m).*dl(n);
                P_prop=P_prop+P*dk(m).*dl(n);
            end
        end
        hrms=(2*P_prop).^.5/2/pi;
        Fr=N.*hrms./ubar;L=min(0.75./Fr,1);
        G1=zeros(5400,2250);G1(~isnan(H2d))=2*L.*L.*E;
        G2d(:,:,a)=G1;
        Fr2=zeros(5400,2250);Fr2(~isnan(H2d))=Fr;
        Fr2d(:,:,a)=Fr2;
        clear G1 G Fr2
        a=a+1;
    end
    fid=fopen(['Eflux_y',num2str(year),'m',num2str(month),'.bin'],'w','b');
    fwrite(fid,G2d,'single')
    fid=fopen(['Fr_y',num2str(year),'m',num2str(month),'.bin'],'w','b');
    fwrite(fid,Fr2d,'single')
end

